A wave-function Monte Carlo method for simulating conditional master equations 
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Wave-function Monte Carlo methods are an important tool for simulating quantum systems, but 
the standard method cannot be used to simulate decoherence in continuously measured systems. 
Here we present a new Monte Carlo method for such systems. This was used to perform the 
simulations of a continuously measured nano-resonator in [Phys. Rev. Lett. 102, 057208 (2009)]. 
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I. INTRODUCTION 

The now standard "wave-function Monte Carlo 
method" for simulating the evolution of a quantum sys- 
tem undergoing decoherence is a very important numer- 
ical tool [JQ. This method allows a simulation of the 
density matrix, an object of size N'^ where N is the di- 
mension of the system, to be replaced by a simulation 
of a number of pure states, each of which is only of size 
N. With the increasing relevance of continuous measure- 
ment and feedback control 0] to experimental quan- 
tum systems, especially in superconducting circuits [^,0 
and nanomechanics [10l4l4l |. one needs to simulate con- 
tinuously measured systems subject to decoherence. The 
standard Monte-Carlo method cannot be used in this 
case, because it applies to master equations but not to 
stochastic (or conditional) master equations (SME's). 

To date, two Monte Carlo methods have been devised 
for simulating conditional master equations, but both 
suffer limitations. The first is by Gambetta and Wise- 
man [l^ , who used the linear formulation [l^, [13] of 
quantum trajectories to derive their method. The less de- 
sirable feature of this method is that it requires evolving 
a fraction of ensemble members that end up contributing 
negligibly to the final density matrix, and to this extent 
it is inefficient. The second method, recently suggested 
by Hush et al. [3, is specifically designed for simulat- 
ing systems with very large state-spaces, in which it is 
not possible to use wave-function methods. This requires 
the use of a quasi-probablity density, such as the Wigncr 
function, and is therefore not as simple to apply to many 
systems. Further, the elements in the ensemble for this 
method are not wave- functions but points in phase space. 
This is important for very large state-spaces, but less de- 
sirable when wave-functions (pure-states) can be used. 
Here we present a wave-function Monte Carlo method 
that avoids all the above issues. This method was used 
to perform the simulations in reference [l9j . but the de- 
tails were not presented there. 

In the next section we state the standard Monte Carlo 
method for reference purposes. In section Hill we present 
the new method with a minimum of discussion. The 
purpose is that this section should serve as an easily ac- 
cessible reference for anyone wanting to implement the 
method. We also note that a parallel implementation us- 



ing C++/MPI is available from the author's website |20| . 
In section IIVI we show how the method is derived, and 
thus show that it reproduces the evolution of a stochas- 
tic master equation. In section |V] we use the method to 
simulate a measurement of the energy of a harmonic oscil- 
lator, and compare it to a direct simulation of the SME. 
Section [VII concludes with a summary of the results. 



II. THE STANDARD MONTE CARLO METHOD 

In what follows. L and A/ are operators, p is the density 
matrix, and dW is a Wiener process, independent of any 
other Wiener processes that may be introduced. 

The standard wave-function Monte Carlo method is 
implemented as follows(see, e.g. [3|). To simulate the 
master equation 



= -7(LtLp + pL^L - 2LpL'i) 



(1) 



we perform the following steps: 

1. Create a set of TV pure states \tpn), so that the 
desired initial value of p is approximately 



1 ^ 



(2) 



2. Evolve each pure state by repeating the following 
steps (i and ii): 

i) Increment each state using the stochastic Schrodiner 
Equation (SSE) 

d\i>n) = -7 [L^ -2{L + L^)J L\i'n)dt 

+ ^L\ijn)dVn, (3) 



where 



(4) 



and the dVn are mutually independent Wiener noise in- 
crements satisfying (dVn)'^ = dt. 

ii) Normalize each of the iV'n)- 

3. The density matrix at time t is (approximately) 



1 ^ 



(5) 
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III. THE NEW MONTE CARLO METHOD 

The conditional (stochastic) master equation 

dp ^ --^{L^Lp + pL^L - 2LpL'^)dt 

-kiM'^Mp + pMH'I - 2MpAP)dt 
+V2k{Mp + pM'f - {M + M'')p)dW (6) 

describes a measurement of M, and decoherence due to 
an interaction with L. 

To simulate the above SME we perform the following 
steps: 

1. Create a set of N pure states, and N probabilities 
Pn , so that the desired initial value of p is approximately 



N 



N 



p(o)==^p„|^„)(V„|, E^" = i- 



(7) 



n = l 



n=l 



Since the P„ are the weightings of the pure states in the 
ensemble that forms p, the effective size of the ensemble 
is no longer N, but can be characterized, for example, by 
the exponential of the von Neumann entropy of the set 
{Pn}: 



iVeff = exp 



N 



< N. 



(8) 



This effective size is maximized (equal to N) iff all the 
Pn are equal to We therefore choose Pn = l/N as 

the initial values of the weightings. 

2. Evolve each pure state by repeating the following 
steps (i - vii): 

i) Increment each state using the SSE 

d|V„) = -7 [L^ -2{L + Lt) J i|V„)dt 

+ ^L\^Pn)dVn, (9) 



where 



(10) 



and the dVn are mutually independent Wiener noise in- 
crements satisfying (dVn)'^ = dt. 

ii) Normalize each of the \4'n)- 

iii) Increment each state by 



+ ^M\'^Jn)dW, 



where 



N 



{M + A/t) = J2 PrMn\{M + A/t)|Vj, 



iv) Update the probabilities P„ using 



(11) 



(12) 



v) Normahze the P„: P„ Pn j YZ=iP-n.- 

vi) Normalize each of the IV'n)- 

vii) Every few iterations perform the following opera- 
tion (which might be referred to as "splitting", "breed- 
ing", or "regenerating" the ensemble): For each pure 
state whose probability Pj is less than a fixed thresh- 
old Pthrcsh <C 1, we pick the state from the ensemble, 
\ipm), whose probability, P^, is currently the largest in 
the ensemble. We then set j^/'j) equal to \ijjm), thus eras- 
ing |-0j) from the ensemble. We set both P, and Pm equal 
to Pml'2- Thus we have "split" the highest probability 
state into two members of the ensemble, and this state is 
(most likely) no longer the highest contributing member. 
After we have done this for each Pj < Pthrcsh, we then 
normalize all the P„ as per v) above. 

3. The density matrix at time t is (approximately) 



N 



(14) 



n=l 



A. Considerations for Numerical Accuracy 

In the standard Monte Carlo method the only param- 
eter that we must chose to reach a desired accuracy is 
A^; we merely increase N until we obtain this accuracy. 
For the new Monte Carlo method we have two param- 
eters that affect the error. The first is the minimum 
effective ensemble size during the evolution, min(iVeff). 
The second comes from the regeneration step. In each 
regeneration we eliminate some states. If we denote sum 
of the probabilities for these "dropped" states as Pdrop, 
then the maximum value of Pdrop during the simulation 
bounds the error from the regeneration step. So to ensure 
numerical accuracy we require that 



min(7Veff) > 1, 
max(Pdrop) < 1- 



(15) 



The values of these two quantities are determined jointly 
by N and Pthrcsh- For a given value of N, there is some 
optimal value of Pthrcsh that ensures that min(A^off) is 
large while keeping max(Pdrop) small. 

For a given simulation it is simple to check whether 
and Pthresh give sufficient accuracy. One merely runs 
the simulation a second time with the same realization 
for the measurement noise dW, and different set of real- 
izations for the noises that model the decoherence, dVi. 
The difference between the two simulations gives one an 
estimate of the error. 



B. Multiple Decoherence Channels and Multiple 
Measurements 



Pn Pn(V'nlV'n). 



(13) 



For simplicity we presented the Monte Carlo method 
for an SME with only a single source of decoherence and 



3 



single measurement. Extending this to m sources of de- 
coherence and I measurements is very simple. A system 
subjected to I continuous measurements and m sources 
of decoherence is described by the SME 

m 

dp^'Y. 7(^1 + PLIL, - 2L,pLl)dt (16) 

i=l 
I 

- kj{M]Mjp + pM]Mj - 2MjpM])dt 

3=1 

I 

+ Y V^iMjP + pM] ~ {Mj + M])p)dW,, 

where the dWj are mutually independent Wiener pro- 
cesses. To simulate this SME one simply repeats steps 
2. i and ii for each of the m decoherence channels, and 
steps 2. iii - vi for each of the / measurement channels. 



C. Inefficient Measurements 

The form of the SME given in Eq. lfTB)) above is general 
enough to include inefficient measurements Q . To make 
the J*'' measurement inefficient we simply choose one of 
the Li to be equal to Mj, and adjust the values of the 
corresponding 7^ and kj to obtain the desired efficiency. 



D. Using Milstien's Method for Time-Stepping 

If one simulates a stochastic differential equation 
(SDE) simply by replacing dt with a small time-step Ai, 
and dW by a zero mean Gaussian random variable with 
variance At (which we will call AVF), then the solution 
is only guaranteed to be accurate to half-order in At. 
Ensuring that the simulation is accurate to first-order 
in At is simple, and the method for doing this is called 
Miltstien's method. Milstein's method involves adding a 
term to the differential equation that is proportional to 
{AW^ - At). The exact form of the Milstien term de- 
pends on the form of the stochastic term in the SDE. If 
the stochastic term is simply a linear operation, the this 
term is given by applying the linear operation twice, and 
multiplying by one half [2l|. Thus, the Milstien term for 
an SDE with the stochastic term aX\ip), for a number a 
and operator X, is 



A|7A)Mil = y(AT^2_^^)^2|^^^ 



IV. DERIVING THE METHOD 



(17) 



We begin by noting that if we apply the part of the 
evolution containing L first, and that containing M sec- 
ond, we get the evolution correct to first-order in dt. If 
the density matrix is given by p = J^n -fn|V'n)j then the 



L part of the evolution is obtained by using the standard 
Monte Carlo method (steps 2. i and ii above). To simu- 
late the part containing M (the measurement part), we 
note that this evolution can be written as lal 



1 



pit + dt) = j^A{a)p{t)AHa), 



(18) 



where A is an operator that depends on the measurement 
result, a, and Af is simply an overall normalization factor. 
The measurement result a is the real number 



a = 2k{M + M^)dt + y^dW. 

By substituting p ^ Y.n -Pri|V'n)(V'«| into Eq.((T 
that 

1 ^ 

p{t +dt) = — Y PnA\^n) {'<Pn\A^ 



(19) 



we find 



U 

N 

E 

n=l 



n=l 



P^{lPn\A^ A\^n) 

Af 



A|^„)(7A„|At 

{^n\A^A\i'n) 



Y Pn{t + dt)\Mt + dt)){lPn{t + dt)l 



n=l 



which gives us the simple update rules 

p^^,^,^^E^mJl^^l^^ (20) 



\lpn{t + dt)) 



A\ip„ 



(21) 



where N is chosen so that ^„Pn{t + dt) — 1. From 
Eq.(29) in reference Q, the operator A is 

A{a) = 1-7 [M^ - 2{M + aP)] Mdt + ^J¥fMdW, 

and this gives us the evolution sequence for the Monte 
Carlo method presented above. 

The effect of the measurement is to increase the prob- 
abilities of some states, and reduce those of others. This 
reduces the effective size of the ensemble, and before too 
long there will only be one state left in the ensemble. To 
correctly model the noise being introduced into the sys- 
tem by the decoherence (the part of the evolution con- 
taining L) we need to maintain a large number of states 
in the ensemble. We solve this problem by using the "re- 
generation" procedure (step 2. vii). Once every so-often 
we discard those states from the ensemble whose prob- 
abilities, and thus contribution, has become negligible. 
This discarding process does not effect the density ma- 
trix unduly so long as the total amount of probability 
of the discarded states is very small. Once the "small" 
states have been discarded, we must choose new states 
to replace them, and we must do this without affecting 
the density matrix. This is easily achieved by duplicat- 
ing some of the states that have a large contribution, and 
dividing the probability for each of these states equally 
between the original state and its duplicate. One should 



duplicate the states with highest probabihty, as this pro- 
vides the biggest increase in the effective size of the en- 
semble. With the addition of this regeneration procedure, 
our Monte Carlo method is complete. 



V. EXAMPLE OF A NUMERICAL 
SIMULATION 

Here we simulate a continuous measurement of the en- 
ergy of a harmonic oscillator, using both the SME and 
the Monte Carlo method. We subject the osciUator to a 
randomly fluctuating white noise force, which serves as 
a simple model of (infinite temperature) thermal noise. 
The evolution due to the fluctuating force is given by [l^l 



(22) 



where /? determines the strength of the force. The term 
in this equation proportional to dt is due to the transfor- 
mation from Stratonovich to Ito noise. Since the observer 
does not know the fluctuating value of the force, she must 
average over it. The evolution of the observer's density 
matrix is then given by the master equation 



dp = — /3[a;, [x, p\]dt. 



(23) 



Adding to this the evolution due the continuous measure- 
ment of energy (equivalently a continuous measurement 
of the phonon number, N = a^a), and the Hamiltonian 
evolution, the SME is 



dp = -iLu[N,p\dt- ^[x, [x,p\]dt - k[N, [N, p\]dt 



2k{Np + pN - 2(N)p)dW. 



(24) 



Here w is the frequency of the oscillator, x = {a + a^) 
is the dimensionless position, and k is the measurement 
strength. 

We simulate this equation using the unnormalized ver- 
sion 

dp= -iu}[N,p]dt- ^[x, [x,p]]dt- k[N, [N,p]]dt 



+ {Np + pN){4k{N)dt + V2kdW) 

+k{dW^ - dt){N'^p + pN^ + 2NpN), (25) 

and normalizing p after each time-step. The reason we 
use this unnormalized version, which you will note is the 
same unnormalized version that we use for the wave- 
function in the Monte Carlo method, is that it makes the 
noise term in the SME linear in p, which in turn makes 
the Milstien term simpler to calculate. In the above equa- 
tion the Milstien term is the final term. Since this equa- 
tion is only accurate to first-order in the noise part of the 
evolution, we have also only included a first-order term 
for the deterministic evolution due to the Hamiltonian. 




FIG. 1. (Color online) The average value of the phonon num- 
ber, A'', for a continuous measurement of A'^, for an oscillator 
driven by a white-noise force, (a) A direct simulation of the 
SME (black line) and a simulation using the Monte Carlo 
method (grey line (blue online)). In this case the time-step 
is df = 2 X 10^*T, and the ensemble has f024 members, (b) 
A zoomed-in version of the Monte Carlo simulation in (a) 
(medium grey line (blue online)); a direct simulation of the 
SME with half the time-step (black line); and the MC simu- 
lation with half the time-step (light grey line). 



The corresponding evolution for the Monte Carlo wave 
function is 



iujN - ^x^ - kN^ 



dt\i:) 



(26) 



-{Ak{N)dt + V2kdW)N\i:) + iy//3xdV[ 



+k{dW^ - dt)N^ 



+ ^{dV^^dt)x^\i;), 



where dV is uncorrelated with dW. 

We now simulate the SME, using both Eq. (p5)) and the 
Monte Carlo method, and compare the results. For this 
simulation we set k = g = O.lf (where / = a;/27r = 
1/T), start the oscillator in the Fock state with three 
phonons. We use a Fock-statc basis, and truncate the 
state space at 9 phonons. For our first run we choose 



5 



dt = 2x and run for a time of t = lOT. For 

the Monte Carlo run we choose the ensemble size to be 
Ncns — 1024, and Pthrosh = 0.2/iVcns- This choice results 
in min(7Veff) = 745.2 and max(Pdrop) = 0.003. We plot 
the expectation value of the phonon number for both 
simulations in Fig. [1^. We see that the results agree, but 
the solutions slowly diverge. To determine the source 
of this divergence wc perform to more simulations. For 
the first one we double the size of the ensemble, and for 
the second we halve the time-step. Note that when we 
halve the time-step, we must use a noise realization that 
is consistent with that used for the first run, so that we 
can directly compare the trajectories in both cases [2l| . 

We find that doubling the size of the ensemble has lit- 
tle effect on the result of the Monte Carlo simulation. 
Halving the time-step, on the other hand, reduces the 
divergence between the two simulations considerably. In 
Fig. [T|d we plot the direct simulation of the SME using 
the smaller time-step, along with the Monte Carlo sim- 
ulations using both time-steps. This plot is zoomed-in 
version of the trajectory in Fig [T|d. These results show 
us that the ensemble size of 1024 is sufficient for this sim- 
ulation, the inaccuracy being due almost entirely to the 
finite size of the time-step. 



ous observation, while also being subjected to noise and 
decoherence. This method is more efficient than the pre- 
viously available method [l5|. We have also applied it 
to an example system, determining in this case sufficient 
resources to reproduce the SME. 

Note added: Upon writing up this work, we discovered 
that a key element, that of "splitting" the ensemble, had 
been introduced previously by Trivedi and Ceperley for 
Monte Carlo simulations of classical systems. See |23j . 
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VI. CONCLUSION 



We have presented a wave-function Monte Carlo 
method for simulating systems that are under continu- 
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